
diffmean_helper <- function(data, treatment, item, m=NULL, m_idx){

  # m: moderator variable name
  # m_idx: modeartor index (0 or 1 if moderator is binary)

  if (!is.null(m)){
    data <- subset(data, data[,m] == m_idx)
  }

  treat_filter <- data[,treatment] == 1 
  control_filter <- data[,treatment] == 0 
  treat <- subset(data, treat_filter) 
  control <- subset(data, control_filter)
  point <- mean(treat[,item], na.rm=T) - mean(control[,item], na.rm=T)
  n1 <- sum(!is.na(treat[,item]))
  n0 <- sum(!is.na(control[,item]))
  se <- sqrt(var(treat[,item], na.rm=T)/n1 + var(control[,item], na.rm=T)/n0)
  upper <- point + qnorm(0.975)*se
  lower <- point + qnorm(0.025)*se
  return(c(point, se))
}

diffmean_block <- function(data, treatment, item, m=NULL, m_idx){

  data_ls <- split(data, data$block_rand)
  group_est <- as.data.frame(t(sapply(data_ls, diffmean_helper, treatment,
                                      item, m=m, m_idx=m_idx)))
  colnames(group_est) <- c("ate", "se")
  block_n <- as.numeric(table(data$block_rand)) # # of obs in each block
  group_est$n <- block_n
  N <- sum(group_est$n[!is.na(group_est$ate)])

  # weighted sum for point estimate
  point <- sum(group_est$ate * group_est$n/N, na.rm=T)
  # weighted sum for se
  se <- sqrt(sum(group_est$se^2 * (group_est$n/N)^2, na.rm=T))
  
  upper <- point + qnorm(0.975)*se
  lower <- point + qnorm(0.025)*se
  p <- 2*(1-pnorm(abs(point/se)))
  return(c(lower, point, upper, p, se))
}
